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Abstract 

We perform error analyses explaining some previously mysterious phenomena aris- 
ing in numerical computation of the Evans function, in particular (i) the advantage of 
centered coordinates for exterior product and related methods, and (ii) the unexpected 
stability of the (notoriously unstable) continuous orthogonalization method of Drury in 
the context of Evans function applications. The analysis in both cases centers around 
a numerical version of the gap lemma of Gardner-Zumbrun and Kapitula-Sandstede, 
giving uniform error estimates for apparently ill-posed projective boundary-value prob- 
lems with asymptotically constant coefficients, so long as the rate of convergence of 
coefficients is greater than the "badness" of the boundary projections as measured by 
negative spectral gap. In the second case, we use also the simple but apparently previ- 
ously unremarked observation that the Drury method is in fact (neutrally) stable when 
used to approximate an unstable subspace, so that continuous orthogonalization and 
the centered exterior product method are roughly equally well-conditioned as methods 
for Evans function approximation. The latter observation makes possible an extremely 
simple nonlinear boundary-value method for possible use in large-scale systems, extend- 
ing ideas suggested by Sandstede. We suggest also a related linear method based on the 
conjugation lemma of Metivier-Zumbrun, an extension of the gap lemma mentioned 
above. 
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1 Introduction 



Recently, numerical Evans function computations have received a great deal of attention 
as a tool for the stability analysis of standing and traveling wave patterns in one and 
several dimensions; see, e.g., [Ml ^ |B?Zl iBDGl EM iHSZl ItHZTl lmIZ2[ iLPSSl iGLZl 
IHLZl IUHNZI HLyZl, HLyZ2| . In the work of the author together with Brin, Humpherys, 



Sandstede, and others, there has emerged a small list of three computational rules of thumb, 
without which Evans function computations become hopelessly inefficient, but with which 
they become in usual situations almost trivial. Indeed, Humpherys has developed a general 
MATLAB-based package (STABLAB) based on these principles that gives excellent results 
on essentially all problems up to now considered. 

Two of the items on this list are self-evident, but the third, the need to "center" coor- 
dinates, does not appear to be well-known and, indeed, at first sight appears to contradict 
standard stability principles. The purpose of this paper is twofold: first, to share these 
practical rules of thumb and, second, to give a mathematical justification for the better- 
than-expected observed results of their implementation, that is, to put these ad hoc princi- 
ples on a rational and quantitative basis. In the process, we discover a numerical analog of 
the gap lemma of [GZl IKSj . a sort of super convergence principle; a new stability property 
of the well-known continuous orthogonalization method of Drury; and, building on ideas 
of Sandstede [S] and Humpherys-Zumbrun [HuZlj . an extremely simple boundary- value 
method for possible use in ultra-large scale systems. 

1.1 Computation of the Evans function 

Let L be a linear differential operator with asymptotically constant coefficients along some 
preferred spatial direction x, and suppose that the eigenvalue equation 

(1.1) {L-X)w = 

may be expressed as a first-order ODE in an appropriate phase space: 

(1.2) W^ = A{x,X)W, lim A{x,X)=A±{X), 

^ ' X— >±oo 

with A analytic in A as a function from C to C^(M, C"^") and the dimension k of the stable 
subspace 5+ of A^ and dimension n — k of the unstable subspace U- of A^ summing to 
the dimension n of the entire phase space. Then, the Evans function is defined as 

(1.3) Z)(A):=det(l^+ ••• W+ W,-^, ••• W-)^^^^, 

where , . . . , and VF^^^, . . . , W^^ are analytically-chosen (in A) bases of the manifolds 
of solutions decaying as x ^ +oo and — oo. For details of this construction, see, e.g., 
[SGJl IFWl KSl EH mi Imlzn IHSZ] and references therein. 

Analogous to the characteristic polynomial for a finite-dimensional operator, D{-) is 
analytic in A with zeroes corresponding in both location and multiplicity to the eigenvalues 
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of the linear operator L |GJlt[GJ2] . Taking the winding number around a contour T = dA C 
{9?A > 0}, where A is a set outside which eigenvalues may be excluded by other methods (e.g. 
energy estimates or asymptotic ODE theory), counts the number of unstable eigenvalues in 
A of the linearized operator about the wave, with zero winding number corresponding to 
stability. See, e.g., [BH [B^ IBDUI IHSZI ImlZTl IBHKZl iHLZl ICHJNZl |HLyZlHHLyZ2t [BHZ] . 
Alternatively, one may use Mueller's method or any number of root-finding methods for 
analytic functions to locate individual roots; see, e.g., |OZl ILS] . 

Numerical approximation of the Evans function breaks into two steps: (i) the compu- 
tation of analytic bases for stable (resp. unstable) subspaces of (resp. A_) and (ii) 
the propagation of these bases by ODE (jl.2p on a sufficiently large interval x £ [M, 0] 
(resp. X € [— M, 0]). In both steps, it is important to preserve the fundamental property 
of analyticity in A, which is extremely useful in computing roots by winding number or 
other methods. Both problems concern (different aspects of) numerical propagation of sub- 
spaces, the first in A and the second in x, thus tying into large bodies of theory in both 
numerical linear algebra |ACR[ IDDFt IDElt IDE21 IDF] and hydrodynamic stability theory 
[raimiNRniNR2l[NE^[NEil [B]. 

Problem (i) has been examined in [HSZl IZ2[ IBHZj : for completeness, we gather the 
(existing but dispersed) conclusions here in Section [6l Our main emphasis, however, is on 
problem (ii) and numerical stability analysis, for which we obtain substantially new results. 

1.2 Three Bad Things: numerical pitfalls and their resolutions 

We now focus on problem (ii). Our three basic principles for efficient numerical integration 
of (|1.2p are readily motivated by consideration of the simpler constant-coefficient case 

(1.4) = AW, A = constant. 

1.2.1 Potential pitfalls 

We note the following three basic pitfalls, two obvious and one perhaps less so. 

1. Wrong direction of integration. Consider the simplest case that the dimension of 
the stable subspace of A is one, so that we seek to resolve a single decaying eigenmode, with 
all other modes exponentially growing with increasing x. Evidently, the correct direction of 
integration is the backwards direction, from x = +M back to x = 0, in which the desired 
mode is exponentially growing, and errors in other modes exponentially decay. Integrating 
in the forward direction would be numerically disastrous, with exponential error growth 
e^*^, where rj is the spectral gap between decaying and growing modes (recall, M is large): 
the analog for Evans function computations of integrating a backward heat equation. In 
general, we must always integrate from infinity toward zero. 

2. Parasitic modes (related to 1). For general systems of equations, the dimension 
of the stable subspace of A typically involves two or more eigenmodes, with distinct decay 
rates Hi < H2 < 0. Integrating in backward direction as prescribed in part 1 above, we 
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resolve the fastest decaying ni mode without difficulty. However, in trying to resolve the 
slower decaying fi2 mode, we experience the problem that errors in the direction of the /ii 
mode grow exponentially relative to the desired /i2 mode, at relative rate e^^^~^^^'^^ . That 
is, parasitic faster-decaying modes will tend to take over slower-decaying modes, preventing 
their resolution. 

Remark 1.1. Degradation of results from parasitic modes is of the same rough order as 
that resulting from integrating in the wrong spatial direction, differing only in the fact that 
the maximum spectral gap between two decaying modes is typically one-half or less of the 
maximum gap between decay and growing modes. Thus, it is crucial to address this issue. 

3. Nonequilibrium state. A more subtle problem is that integrating even a single 
scalar equation w' = —aw, a > 0, over a long interval [M, 0], leads to (sometimes quite large) 
accumulation of errors, and, more important, a large number of mesh points/computations. 
The single exception is the equilibrium case a = 0, which for essentially all numerical ODE 
schemes is resolved exactly. 

This can be understood more quantitatively by the following heuristic computation, 
assuming a perfectly adaptive scheme and no machine error. The truncation error Tj for a 
kth order scheme at step j is proportional to the {k + l)th derivative of the solution times 
Axj , where Axj is the size of the jth step, or Cfca'^"'"^e~"^^ Ax^ . Taking tj ~ TOL for some 
fixed tolerance TOL, we thus obtain Cka''~^^e~'^^^ Axj ~ TOL, or 

Aj ^' 

Inverting, and integrating « ^ from to M0 we obtain an estimate 

f-M 

(1.5) J ~ cl^'^TOL'^/'^ka^/'' / (0/^)6-^^^/^ ~ c^/^OL-i/'^ifcaV^ 

Jo 

as Af — 5- +00 for the total number J of mesh blocks, which goes to zero as a — s- and to 
infinity as a ^ 00. 

Though it is tempting to think of this as an example of numerical stiffness, that is not 
the case, since the problem involves but a single mode. It seems rather to be a secondary, 
previously unremarked, phenomenon, that in usual circumstances would be neglible. For 
Evans function computations, however, we observe a difference in computational efficiency 
of an order of magnitude or more between the cases a = 1 and a = 0, for TOL ~ 10^^. 

Remark 1.2. It would be iteresting to compare (jl.Sp to results for the standard RK^5 
scheme with which most Evans computations have been done applied to the constant coeffi- 
cient scalar problem w' = —aw, in particular the a)-!^ rate. For variable- coefficient systems, 
additional effects having to do with "conjugation errors" appear as well; see Section\^ 



^ Direction of integration is symmetric for a single mode. 
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1.2.2 Solution one: the centered exterior product method 

Problem 1 is easily avoided by integrating in the correct direction. Problem 2 may be 
overcome by working in the exterior product space A • • ■/\Wil G C^*) (resp. W/T^^ A • • • A 

W~ G C^"-fe)), for which the desired subspace appears as a single, maximally stable (resp. 
unstable) mode, the Evans determinant then being recovered through the isomorphism 

(1.6) det(l^+ ••• W+ W,^, ■■■ W-)r.iW+A---AW+)A{W,^,A---AW~); 

see [ASl [13 EEZl iBDGl IAIB] and ancestors [GBl iNRll [NE2l INRHI [NR4] . This reduces the 
problem to the case k = 1. Problem 3 can then be avoided by factoring out the expected 
asymptotic decay rate e^^ of the single decaying mode and solving the "centered" equation 

(1.7) Z' = {A- fiI)Z, Z(+oo) = r : A+r = 

for Z := e~^^W , which is now asymptotically an equilibrium as x — > +oo. With these 
preparations, one obtains excellent results (BDGl IHuZlj : however, omitting any one of 
them leads to a loss of efficiency of at least an order of magnitude in our experience |HuZ2] . 

1.2.3 Solution two: the polar coordinate method 

Unfortunately, the dimension (^) of the phase space for the exterior product grows ex- 
ponentially with n, since k is ~ n/2 in typical applications. This limits its usefulness to 
n < 10 or so, whereas the Evans system arising in compressible MHD is size n = 15, k = 7 
|BHZ] . giving a phase space of size (2) = 6,435: clearly impractical. A more compact, 
but nonlinear, alternative is the polar coordinate method of [HuZlj . in which the exterior 
products of the columns of W± are represented in "polar coordinates" (J^,7)±, where the 
columns of n+ G C"^'' and G C^""'')^'' are orthonormal bases of the subspaces spanned 
by the columns of W+ := {W^ ■ ■ ■ W^) and W- := {Wj;^^ ■ ■ ■ W') , defined as 
in (jl.3p . i.e., W-^. = r2+a+, W- = f]_a_, and 7± := deta±, so that 

A---A W+A = 7+(17^ A • • • A n^), 

where denotes the jth column of r2±, and likewise WjTj^-^A- ■ - AW' = ^-{Q'La- ■ ■AQ^~''). 
This yields the block-triangular system 

= {i- ^}n*)An, 

^^'^^ (log 7)' = trace {n*An) - trace {n*An){±oo), 

7 := ^e~^^^^'^ {n* An)(±cx))x ^ which the "angular" 0-equation is exactly the continuous 
orthogonalization method of Drury [Drl IDa] . and the "radial" 7-equation, given J7, may 
be solved by simple quadrature. Ignoring the numerically trivial radial equation, we see 
that problem 2 by fiat does not occur. Likewise, for constant A, it is easily verified that 
invariant subspaces Q oi A are equilibria of the flow, so problem 3 does not occur. Indeed, 
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for A constant, solutions of the 7-equation are constant, so that any first-order or higher 
numerical scheme resolves log 7 exactly; thus, the 7-equation may for simplicity be solved 
together with the O-equation, with no need for a final quadrature sweepH The Evans 
function is recovered, finally, through the relation 

(1.9) D{X) = det{W+ ■■■ W+ W^^^ ••• VF„-)U=o = 7+7-det(l^+,J^-)U=o. 
1.3 Further questions and description of results 

The discussion of the previous subsection leaves open some important questions. First, can 
we justify this heuristic discussion with rigorous, quantitative error bounds for the actual, 
variable coefficient problems that occur in practice? In particular, we note that centering 
equations as in Section 11.2.21 goes counter to the intuition afforded by standard two-point 
boundary- value theory on intervals [0, M] as M — > 00 {Bel] , which asserts convergence 
error of order e~^^ where r/ is the minimum spectral gap of decaying (resp. growing) 
modes from zero to the solution on [0, -|-c«). Applied blindly to the centered equations, this 
would predict nonconvergence rather than the good behavior observed in practice. 

Second, there is a well-known problem of instability of the continuous orthogonaliza- 
tion method with respect to perturbations disturbing the assumed orthonormal structure 
of Vl |Dal IBrRej . In the language of [BrRej . the Stiefel manifold 5 := {fi : ^1*0, = I^} of 
orthonormal matrices is preserved by the flow of (jl.Sp . but is typically neither attracting 
nor repelling. In view of Remark 11.11 this should lead to terrible results for Evans func- 
tion computations. This issue was discussed at length in |HuZl| . with numerous different 
solutions discussed, from artificial stabilization to geometric integration. Yet, surprisingly, 
the method that performed best was the original Drury algorithm with no stabilization, 
implemented by a standard RK45 scheme. This yielded results quite similar to those of the 
exterior product scheme, which seems to contradict the conclusions of Remark II. 1[ 

Result I. Our first main result is to establish a numerical version of the gap lemma of 
[GZl IKSj . which states that, ignoring machine error, provided the coefficient matrix A{x, A) 
is uniformly exponentially convergent as x — > -|-oo, with rate \A — A+l < Ce~^^ for x > 0, 
and provided the gap between fi minimum real part of the eigenvalues of A^ is strictly 
greater than —6, then the solution of problem (II. 7p on [0, M] initialized as Z(M) = r 
converges as M ^ cxd to the solution on [0, +00) at rate C{9,6)e~^^^ for any < 9 < 9. 
This resolves the first issue, explaining the observed convergence of the centered exterior 
product method. 

Completing the analogy to |GZj . we establish a corresponding result for general centered 
two-point boundary problems (11.70 with projective boundary conditions on [0, Af], under 
the assumption that the gap 7 between /i minimum real part of the eigenvalues of A^ asso- 
ciated with the projective boundary condition at M is strictly greater than —9, obtaining 
convergence at the same rate C{9,9)e~^^^ as M — > +00 for any < 9 < 9. This could be 
viewed as a type of super convergence, as the standard theory [Bel] predicts convergence at 

^See Section 16.3.21 for numerical prescriptions (n,7)± of (^,7) at ±00. 
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rate e "^^j with a nonpositive spectral gap 7 < corresponding to ill-conditioned boundary 
conditions. For detailed statements and proofs, see Section [3l 

Result II. Our second main result is to explain the apparent contradiction between 
observed good results for the polar coordinate method [HuZl] and the well-known insta- 
bility of continuous orthogonalization [Da\ IBrRej . The simple resolution is that, though 
continuous orthogonalization is in general unstable, it is in the present context stable! 

Heuristically, this is quite simple to see. Intuitively, it is is clear that the stable manifold 
of is asymptotically attracting in backward x under the flow of (jl.Sp for confined to 
the Stiefel manifold 5 = {il : = 1^}, and in fact this is well known (see Section [4.ip . 
Thus, we need only verify that the Stiefel manifold, likewise, is attracting in backward x. 
Defining the Stiefel error £{Q) := — I, we obtain after a brief computation the error 
equation 

(1.10) £' = -£{n*An) - {n*Any£ 

of [HuZlj . Linearizing about the exact solution — > Q-^- , £ ^ and replacing coefficients by 
their asymptotic limits, we obtain a linear equation £' = A+£, where A£ := —£A^ — {A-^^)*£ 
is a Sylvester operator, A^ := {Q*^_A^Q-^-) , with eigenvalues and eigenmatrices aj + a"^, rjr^, 
where Uj and rj are eigenvalues and eigenvectors of A^ . Noting that the eigenvalues of A^ 
are exactly the eigenvalues of ^4+ restricted to its stable subspace, i.e., the stable eigenvalues 
of A^, we find that has positive real part eigenvalues, and so £ decays in backward x. 

That is, the Stiefel manifold is attracting under the backward flow of continuous orthog- 
onalization (repelling under the forward flow) if ^ is a stable subspace of A, an observation 
that previously seems to have gone unremarked. This confirms that, as suggested by nu- 
merical results of |HuZlj . continuous orthogonalization is roughly equally well-conditioned 
as the centered exterior product method. For details and further discussion, see Section |H 

Remark 1.3. An important consequence is that geometric integrators like those suggested 
in IBrRej for general Orr-Sommerfeld applications are probably not worth the trouble for 
Evans function computations, since the Drury method is stable and much simpler to code. 

Remark 1.4. The generalized inverse method, or Davey method Wal^ Q' = {I — ^}Q^)AQ, 
where il^ := (i7*0)~^ri* denotes the generalized inverse, exhibits neutral error growth £' = 0, 
so is often used as a stabilization of the basic continuous orthogonalization method l\1.8\i (i) 
of Drury. In the context of Evans function computations, the behavior is slightly worse, 
however JHuZlf . This can now be understood from the fact that the Drury method actively 
damps errors in the Evans function context. The fact that the Drury method outperformed 
the Davey method (and all others) was a mysterious aspect left unresolved in IHuZlf . 

Result III. For equations arising in complicated physical systems or through transverse 
discretization of a multi-dimensional problem on a cylindrical domain [LPSS] , the dimension 
n can be very large. For example, for the multidimensional systems considered in [LPSS], 
n = 8M ~ 48 (see p. 1447, [LPSS] ) where M ~ 6 is the number of transverse Fourier 
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modes being computed. The development of numerical methods suitable for efficient Evans 
function computations for large systems has been cited by Jones and others as one of the 
key problems facing the traveling- wave community in the next generation p]. 

For large dimensions, the accumulation of errors associated with shooting methods ap- 
pears potentially problematic, and so various other options have been considered. For 
example, one may always abandon the Evans function formulation and go back to direct 
discretization/Galerkin techniques, hoping to optimize perhaps by multi-pole type expan- 
sions on a problem-specific basis. However, this ignores the useful structure, and associated 
dimensionality reduction, encoded by existence of the Evans function! 

Alternatively, Sandstede jS] has suggested to work within the Evans function formu- 
lation, but, in place of the high-dimensional shooting methods described above, to recast 
(jl.2p as a boundary-value problem with appropriate projective boundary conditions, which 
may be solved in the original space C"" for individual modes by robust and highly-accurate 
boundary- value/continuation techniques. 

Problems with this scheme as conventionally implemented in uncentered coordinates are 
two. First, solutions exponentially decay as x ^ ±00, so that the direct connection to data 
at 00 of (jl.7p is lost; as a consequence, up to now, it is not known how to recover analyticity 
of the Evans function by such a scheme. Second, since the uncentered problem involves 
decaying modes as well as growing modes, there must be provided boundary conditions 
at X = as well as at a; = itM; indeed, it is the boundary conditions at x = that 
mainly determine the decaying modes we seek. Since behavior of (jl.2p is only known near 
its asymptotic limits as x — s- ±00, there appears to be no analytic way to prescribe a 
priori well-conditioned projective boundary conditions at x = (or, as mentioned already, 
to relate these to a desired asymptotic behavior as x ^ ioo), and so apparently these 
must be adjusted "on the fly" by trial and error, a process that requires error checks and 
additional complications in program structure. 

As pointed out in [HuZl j (last sentence of introduction), using the polar coordinate 
method- or any centered scheme- as the basis of a boundary- value scheme eliminates im- 
mediately the first problem, of preserving analyticity, since in the centered format date is 
explicitly described as x — > ±00. The second problem in general remains. However, by the 
remarkable stability property recorded in result II, the polar coordinate method involves 
only modes that are neutral or growing as x ^ ±00, and not decaying; that is, it is both 
centered and one-sided. The centered exterior product method though dimensionally unsuit- 
able shares this property as well; indeed, it is precisely the one-sided property that makes 
these schemes suitable for shooting. For such schemes, appropriate projective boundary 
conditions by result I are full Dirichlet conditions as x — > ±00, with no boundary conditions 
at X = 0, and thus the second, essentially logistic, problem does not either arise. 

We therefore propose the polar coordinate method with Dirichlet conditions at x = itM 
as a promising candidate for boundary- value-based Evans function computations, noting in 
particular that it is essentially trivial to program given an existing shooting code. The only 
disadvantage that we immediately see is the nonlinearity of the scheme. We propose at the 

See however [GLZj . which points out a useful intermediate structure based on Fredholm determinants. 
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same time an alternative linear, centered but two-sided, scheme based on the conjugation 
lemma of [MeZl [ZT] , an extension of the gap lemma that is our main tool in the analysis. 



1.4 Discussion and open problems 

We gather in this paper a complete prescription together with rigorous error bounds for 
efficient numerical Evans function computations by the shooting methods of |HuZH IHLZ[ 
HLyZl, iBHZj . etc., of systems up to the intermediate size n ~ 20 or so encountered in 



one- and multi-dimensional problems of continuum mechanics, and propose some promis- 
ing boundary- value methods for further exploration in computations for ultra-large systems 
of size n ~ 50 and up. In the process, we rehabilitate the continuous orthogonalization 
method of Drury, explaining its unexpectedly good performance in the context of Evans 
function computations. Finally, and most important, we point out the importance of cen- 
tered coordinates, in contrast to standard numerical intuition and practice in the study 
of boundary- value problems on unbounded domains [Bel], establishing the related stabil- 
ity/super convergence principle embodied by our numerical gap lemma. 

The latter result seems to suggest larger implications in the construction of numeri- 
cal boundary- value schemes. At the same time, it serves to clarify some up-to-now rather 
confusing existing results. For example, in the seminal work [Er], Erpenbeck performed 
a numerical Evans function analysis of stability of ZND detonation waves by a method 
obeying principle 1 and 2 of Section 11.2. 1[ Much later, Lee and Steward [LSj introduced 
what is now the effective standard method in detonation literature, obeying principle 2 
but violating principle 1, reporting an apparently counter- intuitive improvement in speed 
of computation. The explanation of this paradox is that Erpenbeck carried out his com- 
putations in uncentered coordinates, whereas Lee and Steward, by mapping to a bounded 
interval and applying a singular integral solver effectively centered their equations, factoring 
out the principal dynamics. Thus, there is a cancellation of errors involved, that cannot be 
seen without reference to principle 3. A centered version of Erpenbeck's original method 
appears to outperform both schemes by an order of magnitude; see |HuZ2j for further 
discussion / simplifications. 

The main mathematical interest of our numerical convergence results is their application 
to two-sided boundary- value schemes posed on the entire interval [0, M]. Indeed, for shoot- 
ing methods, a routine translation into the discrete setting of the continuous gap lemma 
yields the result, whereas the general case involves a more subtle argument based on approx- 
imate conjugation to constant-coefficients; see Remark 13. 121 The determination of realistic 
mesh requirements for centered boundary- value schemes, and the question of whether or 
not centered boundary- value schemes yield in practice the same good performance observed 
for centered shooting schemes, remain important open problems. 
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2 Preliminaries: the gap and conjugation lemmas 

We begin by recalling the standard gap and conjugation lemmas of \GZ\ IKS] and |MeZ] . 
Consider a general family of first-order ODE 

(2.1) W'-A(x,A)W = 

indexed by a parameter A G $7 C C", where W G C^, x € M and ""' denotes d/dx. Assume 
(ho) Coefficient A(-,A), considered as a function from O into C^{x) is analytic in 
A. Moreover, A(-,A) approaches exponentially to limits A± as x — > ±oo, with uniform 
exponential decay estimates 

(2.2) \{d/dx)''{A - A±)| < Ce^^l^l, for x ^ 0, < A; < i^, 
C, 9 > 0, on compact subsets of 

Lemma 2.1 (The gap lemma [GZ, ZH]). Assuming (hO), ifV^{A) is an eigenvector of A_ 
with eigenvalue /i(A), both analytic in A, then there exists a solution of (j2.ip of form 

(2.3) W(A,x) = y(x,A)e'^(^)^ 

where V is in x and locally analytic in A and, for any fixed 6 < 6, satisfies 

(2.4) y(x,A) = y-(A) + 0(e-^>l|y-(A)|), x < 0. 
Proof. Setting W(x) = e^^V{x), we may rewrite W = AW as 

(2.5) V = {A. - fiI)V + QV, e := (A- A_) = O(e-^l^l), 

and seek a solution V{x,A) V^{x) as x ^ oo. Choose 9 < 9i < 6 such that there is 
a spectral gap j5ft((TA_ — (^ + 6*1))! > between (tA_ and /x + ^i. Then, fixing a base 
point Aq, we can define on some neighborhood of Aq to the complementary A_-invariant 
projections P{A) and Q{A) where P projects onto the direct sum of all eigenspaces of A_ 
with eigenvalues /i satisfying K(/i) < 3?(/x) + 0i, and Q projects onto the direct sum of 
the remaining eigenspaces, with eigenvalues satisfying 5R(/i) > 5?(/i) + 9i. By basic matrix 
perturbation theory (eg. [Kat]) it follows that P and Q are analytic in a neighborhood of 
Aq, with 

(2.6) ^iA-^^.I)xp < de^'''), x > 0, e^^—'^^^'^Q < C7(e^i^), x < 0. 
It follows that, for M > sufficiently large, the map T defined by 

TV{x) = V-+ r e(*— '^■^)(^-^)PG(y)y(y)dy 

(2-7) „_M 

- / e^^—^'^^^^-y^QQiy)V{y)dy 

J X 
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is a contraction on L°°(— cxo, —M]. For, applying (j2.6p . we have 

(/•X p—M 

^ ' eix -{e-ei)M i 

= C\Vi - Faloo < ^\Vi - V2\oo. 

(7 — (71 Z 

By iteration, we thus obtain a solution y € L°°{—oo, —M] oiV = TV with V < C3\V~\; 
since T clearly preserves analyticity V{A, x) is analytic in A as the uniform limit of analytic 
iterates (starting with Vq = 0). Differentiation shows that F is a bounded solution of 
V = TV if and only if it is a bounded solution of (12. Sp . Further, taking Vi = y, V2 = in 
(j2.8p . we obtain from the second to last inequality that 

(2.9) \V -V-\ = \T{V) - T(0)| < Cse^'^lFl < cJ''\V~l 

giving (j2.4p . Analyticity, and the bounds (j2.4p . extend to x < by standard analytic 
dependence for the initial value problem at x = —M. □ 

Remark 2.2. The title "gap lemma'' alludes to the fact that we do not make the usual 
assumption of a spectral gap between /i(A) and the remaining eigenvalues o/A_, as in stan- 
dard results on asymptotic behavior of ODE JCof : that is, the lemma asserts that exponential 
decay of A can substitute for a spectral gap. 

Remark 2.1. In the case ^a{A^) > —9, we may take Q = % and 6 = 9, improving (12. 4p . 

Corollary 2.3 (The conjugation lemma |MeZj ). Given (hO), there exist locally to any 
given Aq € invertible linear transformations P^{x,A) = I + 0+(x,A) and P_(x,A) = 
I + 6_(x, A) defined on x > and x < 0, respectively, ^± analytic in A as functions from 
Q to C'^[0,±oo), such that: 

(i) For any fixed 0<9<9andO<k<K + l, j>0, 

(2.10) \{d/dAy{d/dx)''e±\ < C{j, fc)e-^>l for x ^ 0. 
(a) The change of coordinates W =: P±Z, F =: P±G reduces (j2.ip to 

(2.11) Z'-A±Z = G forx^O. 
Remark 2.2. Equivalently, solutions of (|2.ip may be factored as 

(2.12) w=(/ + e±)z±, 

where Z-t satisfy the limiting, constant- coefficient equations ()2.1ip and Q± satisfy (I2.10p . 
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Proof. Substituting W = P_Z into ()2.ip . equating to (j2.1ip . and rearranging, we obtain 
the defining equation 



(2.13) PL = A_P_ - P_A, P^ ^ I as x 



-oo. 



Viewed as a vector equation, this has the form PL = AP-, where A approaches exponen- 
tially as X — > — oo to its limit A-, defined by 

(2.14) ^_P:= A_P-PA_. 

The limiting operator A- evidently has analytic eigenvalue, eigenvector pair fJ- = 0, P- = /, 
whence the result follows by Lemma |2. II for j = k = 0. The x-derivative bounds < A; < 
K + 1 then follow from the ODE and its first K derivatives, and the A-derivative bounds 
from standard interior estimates for analytic functions. Finally, invertibility of P_ follows 
for X large and negative from ()2.10p and for x < by global existence of a solution to 

(p-i)' = = A+P~^ - P'^A. 

A symmetric argument gives the result for P+. 

□ 

3 A numerical gap lemma 

We now establish the main result of the paper, a numerical analog of Lemma |2. II 

3.1 Continuous problem 

Consider similarly as in ()2.ip a first-order ODE 

(3.1) W'-A{x)W = 0, WeC^ 
on the half-line x € [0, -|-oo), assuming exponential convergence 

(3.2) \{d/dx)''{A -A+)\< Ce-^^, for x > 0, < A: < iiT, 
C, ^ > 0, of A to Aj^ and asymptotic stationarity 

(3.3) y+ G Ker^+. 
Suppose further that there holds the following ga'p condition. 

Assumption 3.1. Yj^ is a k-dimensional invariant subspace of containing all eigen- 
modes associated with nonnegative real part eigenvalues, and no eigenmodes with real part 
< —9, with associated eigenprojection 11+ . 

Let Hq be an arbitrary projection of rank (n — k). 
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Lemma 3.1. For generic Hq, specifically those satisfying (|3.10p below, there exists for any 
a E Range Hq a unique solution of (jS.ip under the projective boundary conditions 



(3.4) noVK(O) = a, lim U+{Wix) - V+) = 0, 

x— >+oo 

satisfying for any < 6 < 9 

(3.5) \{d/dx)''{W{x) -V+)\< C{9,9)e-^'', foi x > 0, < k < K. 
Proof. By Lemma 12.31 there exists an invertible coordinate transformation 

(3.6) w = {i + e)z 

with 

(3.7) \{d/dx)''e\ < C(A;)e~^>l for x > 0, 
converting (13. ip . ()3.4p to 

(3.8) Z' - A+Z = for X > 
and 

(3.9) no^(O) = a, lim U+{Z{x) - V+) = 0, 

X—> + QO 

where Uq := (/ + e{0))-'^Uo{I + 6(0)) is again arbitrary and d := (/ + e(0))"ia. By 
inspection, there exists a unique solution if and only if KerLIo is transverse to i.e., there 
holds the generically satisfied Evans/Lopatinski condition 

(3.10) detnoS+/0, 

where S+ is the complementary (rank n — k) invariant subspace of S^., consisting of the sum 
of the constant solution Z = Vj^ and the solution of the Cauchy problem Z{Q) = a — noV+ 
under the flow of the constant-coefficient equation (j3.8p restricted to the subspace S+ of 
exponentially decaying solutions. □ 



3.2 Discrete problem 

We now consider a discretized version of (|3.ip . (|3.4p on the truncated domain x € [0,M] 
with the corresponding projective boundary conditions 

(3.11) IloW{<d) = a, n+(l^(M)-y+) = 0, 

and examine convergence as M — > +00 and mesh size goes to zero of the approximate to 
the exact solution. 

Remark 3.2. In view of the discussion of the introduction, our main interest in the context 
of Evans function computations is the "Dirichlet" case = arising for one-sided 
schemes suitable for shooting methods: in particular, the centered exterior-project method or 
the polar coordinate method linearized about a desired exact solution. In this more favorable 
case, there is no boundary condition at x = 0, and no arbitrary projection Hq, and the 
projective boundary conditions ()3.1ip reduce to the simple Dirichlet condition W{M) = V+. 
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3.2.1 Difference scheme 



We assume a general linear difference scheme of form 



(3.12) 



SW = 



with boundary conditions 



(3.13) 



no>Vo = a, n+(Wj-y+) = 0, 



where W = (Wi, . . . , Wj) are approximations of W at mesh points Xj, j = 0, . . . , J and 
iS is a linear difference operator with finite stencil {j — i,j + i}, i.e., the value of (5VV)j 
depends on W only through {Wj-i, . . . , Wj+f}. Moreover, we assume that {SW)j depends 
on A linearly and only through the restriction of A to the interval [xj-g, xj^i]. We define 
boundary and truncation errors Eq, e j and r = (tq, . . . , rj) as usual by 



where VVj := W{xj) denotes the solution of the continuous problem (j3.ip . (|3.4p on [0, +oo) 
sampled at mesh points Xj, and hj is the jth mesh length, defined as the maximum difference 
between points Xk involved in the evaluation of Sj . 

This could be a boundary-value scheme, with S realized as a large Jn x Jn banded 
matrix. Or, in the case of our main interest S = C" that boundary conditions are imposed 
only at one end x = M, it could be a backward Cauchy solver such as RK45, with S realized 
as a series of J — i successive {2£ + l)n x (2£ + l)n matrix multiplications. 

3.2.2 Basic assumptions 

We do not specifiy the details of the scheme, other than to make certain mild assumptions. 
Assumptions 3.2. (i) k-th order consistency: As mesh size hj — > 0, 



(a) Strict constant- coefficient stability: In the constant- coefficient case A = A^, if 
^a{A^) < fi, then, for any fi > ix, Cauchy information SW = rh and Wo = Sq imply 
(under appropriate mesh restrictions) 



(3.14) 



SW = hjT, noWo-a = eo, U+{Wj - V+) = ej, 



(3.15) 



Tj \ < Ch] sup \dl-^^W{z) \ 0. 



(3.16) 




i—r<m<j 



0<k<j 



for < r < 2, where A is the backward difference operator {Af)j := fj — fj-i- Likewise, i, 
5ft(T(^+) > ^, then, for any fi < fi, SW = rh and Wj = ej imply 
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In the case S = C" of our main interest, we require only (j3.17p and only for fi < Oo 

As is customary, we ignore machine error. 

Remark 3.3. Condition (i) is of course standard. In the situation of our main interest 
of a backward Cauchy solver in the case S = C", condition (ii) (in this case, the single 
condition (I3.17P ) for r = is equivalent by DuhameVs principle/ variation of constants to 
the homogeneous stability condition 

forxj < Xfe, /i < 0, when SW = and Wk = Sk, which is related to a circle of ideas including 
A-stability and one-sided Lipshitz bounds JDA \DS\ \HNWl\ \HNW^ regarding approximation 
of stiff ODE. For general schemes it may impose impractically severe limitations on the mesh 
size; however, it is satisfied without such conditions for RK45 and other "nice" schemes, 
as may be seen by applying a constant linear coordinate transformation A ^ A := SAS~^ 
for which ^A := (1/2) (A + A*) > jl > ^ with fi arbitrarily close to fi, then applying the 
results of JD1\ \DS\ \IINW1\ \HNW^. For some simple examples, see Remark \3.^ Likewise, 
for Cauchy solvers (ii) (r = 1,2) follows using the difference scheme from (ii) (r = 0). For 
general boundary-value schemes, (ii) must be checked on a scheme-by-scheme basis. For 
general theory on stability of boundary-value schemes, see, e.g., JKri \Bel\ \Be^ . 

Remark 3.4. Evidently, (ii) implies also the dual bounds 

(3.18) iw.l <C(eoe^^^ + sup |/i^r' J] e'^^^^-^'^V^/ifc) 

j-r<m<j 

and 

(3.19) \w^\<c{eje^^^^-^-'^+ sup \h^\' ^ ef^^^^-^^^nh^) 

for differenced data SW = IV {tK), < r < 2, where A denotes the forward difference 
operator {Af)j := fj+i - fj. 

3.3 Discrete conjugation error 

We now make the key observation that the conjugating transformation for the continuous 
problem, up to a small commutation error, is a conjugator for the discrete problem as well. 

Example 3.5. Consider the first-order (forward exphcit/backward imphcit) Euler scheme 

{SW)j := Wj+i - Wj - hjAjWj = 0. 

* This is all that is needed in any case for (I3.17|l . but this seems to give no advantage in the general case 
since there remains the more restrictive assumption (|3.16|l in the other direction. 
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Substituting Wj =: Pj2j, Pj := P{xj), where P = I + is the continuous conjugator of 
(13.61) , we obtain after a brief calculation 



{SZ)j := Zj+i - Zj - hjA+jZj = {"ifZ) 



where S is the realization of the same first-order Euler scheme to the constant-coefficient 
case A = and ^ is the commutator error 

i^Z),) := -h,p;' - P, - hMjP, - P,+iA+,))z,. 



Noting that Vj+i-Vj = hj{AjVj-Vj+iA+j) is a discretization of the ODE P' = AP-PA+ 
defining P, (j2.13p . we find similarly as in (|3.15p that the exact solution P has truncation 
error 

fj := hj' - P, - h,{A,Pj - Pj+iA+j)) = 0{\hj\ sup \P"\) = 0{\h,\e-~^^), 

xe[xj_i,Xj+i] 

we find that 

(3.20) < C\hj\^e-~^''^\Zj\. 

Remark 3.6. Note that the righthand side of p.20p is smaller by a factor hj than the 
corresponding estimate 

{SW)j=0{\hj\ sup \Am-A^+\) sup < C|/i,|e-^"^ sup \Wm\. 

j-e<m<j+e j-i<m<j+e j-^<m<j+l 

obtained by separating out constant- coefficient and exponentially decaying parts in the orig- 
inal equation for W. This additional factor is crucial in the contraction argument of §3.^[ 



Example 3.7. Consider the second-order (forward/backward implicit) midpoint scheme 

(3.21) (5W),+i := W.+i - - (x,+i - x,_i)A,W, = 0. 

Substituting Wj =: PjZj, P = / + © as in (13. 6p . we obtain 

{SZ)j := Zj+i - Zj^i - {xj+i - xj^i)A+jZj = {^Z)j + A{^^Z)j 

where S is the constant-coefficient realization of S, A is the forward difference operator 
(A/)i := fj+i - fj, and 

{^Z)j) := -{xj+i - x,_i)Pri(p,.+i - P,._i - {xj+i - x,_i)(^,P,- - P,+i^+,))z,- 
+ Pf' (^±4^ - pA (.,+1 - x,.r)A,Z, 



{^'Z)j := (pri(P,.+i -P,„i))(Z, 
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Arguing as in the previous example, we obtain similarly as in (|3.20p the bound 

(3.22) |(^Z)j| < C|/ijf e"^^'^' sup \Z^\. 

j-l<m<j 

and, likewise, 

(3.23) |(^^Z)j| < C|/ij|e~^^J sup 

j-l<m<j 

Remark 3.8. The implicit backward Euler scheme of Example VJ. 51 is A-stahle as a Cauchy 
solver in the backward x direction, hence satisfies Assumption \3.2\ in the one-sided case 
= C", /U = 0, with no assumption on the mesh size. For example, in the scalar case 

(3.24) W' = -aW, a > 0, 

yVj+i = yVj + hjaWj yields immediately Wj = (/ + hja)^^Wj-\-i < Wj+i, for any hj > 0. 
The implicit Midpoint method viewed as a two-step backward scheme has characteristic roots 
—ahj lb (ahj)"^ + 1 that remain strictly ^ 1 independent of hj > 0, but is not backward 
A-stable as a Cauchy solver. 

Generalizing the results of the examples, we make the following final assumption on the 
scheme, which appears to be satisfied in most if not all cases of interest. 

Assumption 3.3. Under the change of coordinates Wj =: PjZj, P defined as in ()3.6p . 
equation (13. ip becomes 

(3.25) SZ = ^^Z + /^{^^Z), 

where S is the same discretization scheme used for W applied to the constant- coefficient 
case A = Aj^, A is the forward difference operator, and 

(3.26) \{^''Z)j\<C{e,e)\hj\^-'^e-'^''^ sup \Zm\. 

j~i-l<m<j+e 

3.4 Numerical convergence lemma 

With these preparations, it is straightforward to establish our following main result. 

Theorem 3.9. Assuming (jS^j) . (|3T^ . and Assumptions\3J\\EM and\3M for fixed 

< 9 < 9, for supj \hj\ sufficiently small and M > sufficiently large, there exists a unique 
solution W of ()3.12p . ()3.13p . which, moreover, satisfies for C > independent of M , r, 

(3.27) \Wj - W{xj)\ < C{e-^^^^ + sup \Tj\) for all j = 0, . . . , J, 

i 

where W is the solution of (|3.ip . (j3.4p guaranteed by Lemma \3.1\ and 9 is as in ()3.2p . For 
S+ = C", the mesh condition supj \hj\ « 1 can be relaxed to sup^ 
where —9<i)<Qis less than the smallest real part of the eigenvalues of A+. 
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That is, we assert existence and convergence so long as the spectral gap min^a [A^\■£_^_^ 
associated with the boundary projection 11+ at +00 is greater than —9, where 6 as in (j3.2p 
is the exponential rate of convervence of A to as x — > +00. On the other hand, standard 
theory [Belj requires a positive spectral gap /3 and concludes convergence at rate Ce~^^ 
for Q < (3 < (5. In other words, though it fails the positive gap condition of standard theory, 
the problem remains numerically well-conditioned so long as "badness" of the boundary 
condition as measured by negativity of the spectral gap is less than the rate of exponential 
convergence, in exact analogy with the gap lemma of continuous theory [GZl IKSt IZHj . 

Proof. By Assumption [331_ 5^ = Y.I=q A^^''^ and SZ = J2l=o ^^^""Z^T, where =: 
PjZj, VVj = W{xj) =: PjZj, and hjTj =: Tj. Defining the convergence error £ := Z — Z, 
we thus obtain the error equation 



2 



(3.28) S£ = y"^ /:^''^''£ - 

with boundary conditions 



r=0 



-dM\ 



(3.29) liof'o = and n+(£'j) = -U+{Zj - V+) = 0{e 

Denoting by ^ = r{£) the solution of S£ = Y.l=o ^'''^''£ - ^, dS^S]), we thus obtain 
the fixed-point formulation 

(3.30) £ = T{£) 

for the solution of (j3.28p - (j3.29p . and thus for the solution Z = Z+£ of the original problem. 

Applying bounds (13^8)1 - (|XT9l) of Remark [331 together with ^M), dS^S]), and (I3T5D . 
and the principle of linear superposition, we obtain 

l^il <C(eoe''"^ + J] e^(^^-^'=)(|(M/Of)fe|+ sup \hmm'£)k\) 

0<k<j J-l<m<j 



(3.31) 

+ C(eje'^(-^--^)+ ^ e'^(^^-^^)(|(vl/0f)fc|+ sup \hmm^£)k\) 

j<k<J i-l<m<j 

j<k<J 

vj/ := (^o^^i^^2^^ where /3 < is greater than the largest real part of the eigenvalues of 
associated with the invariant subspace S+ complementary to and —6 < D < v < 
is less that than the smallest real part of the eigenvalues of A^ associated with S4.; see 
Assumption 13.11 From (I3.3ip , we readily obtain by a discrete version of calculation (12. 8p in 
the argument of Lemma 12.11 the a priori estimate 

(3.32) sup < C(e~^*''^ -|- sup ItjI -|- sup sup \£j\), 

o<j<J 3 3 0<j<J 
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which for supj \hj\ sufficiently small implies 

(3.33) sup \£j\ < C{e^^^' + sup |r, |), 

o<i<J j 

hence, by boundedness of \Pj\ and Wj — W{xj) := PjSj, the desired bound (j3.27p . A similar 
estimate yields contractivity of T in the sup-norm, hence existence and uniqueness in the 
class i°°[0,J]. 

Finally, observing in case $]_|_ = C" that only the Y2j<k<j terms in (j3.3ip appear, we 
may bound supj_i<„<j |/im||(^'^^)fc| using j < k by 

sup \hm\ sup e~^''"^\hm\< sup (e(~^~'^)^^|/ij|) sup e'^^"^|/im| 

j-l<m<j k-e~l<m<k+l 0<j<J k~l-l<m<k+t 

and similarly for to obtain by the same argument 

sup < C7(e-^'^-^ + sup|Tj| + sup {e'^-'^~'^^''^\hj\) sup \Ej\), 

0<j<J j 0<j<J 0<j<J 

in place of (j3.32p . yielding existence, uniqueness, and convergence under the relaxed mesh 
condition sup^- \hj\e^~^~'^^^3 << 1. □ 

Remark 3.10. By Remark \3.8l Theorem \3.9\ applies in particular to the one-sided, or 
Dirichlet, case = C", with the implicit backward Euler method of Example 1 5'. 51 

Remark 3.11. Since —9 — i) < by Assumption \3.1\ in the case S+ = C" of our main 
interest, the mesh and truncation error conditions 

sup l/ijle^^""''^''-' << 1 and sup|Tj | < sup l/ijl^e"^""^ << 1 
j j j 

required for convergence allow for arbitrarily large step size as xj +oo, hence our result 
indeed applies to shooting-type schemes with adaptive step size such as are used in practice. 

Remark 3.12. The method of proof in the argument of Theorem \3.4l based on Lemma \2.3l 
is designed to handle the general boundary-value case. For shooting methods, there is a much 
simpler proof on [L,M] with 1 « L < M, as in the proof of Lemma \2.1l just separating 
constant- coefficient and exponentially decaying parts of A in the original W equation and 
obtaining contraction as in the proof of Lemma \2.1\ from largeness of L. The extension to the 
full domain [0, Af] then follows similarly as in the continuous case by standard convergence 
results for the Cauchy problem on a bounded domain. The conjugation argument is used to 
obtain contraction for globally defined schemes on the full domain [0,M]. 
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3.5 Convergence of the centered exterior product method 

The centered exterior product method as described in the introduction consists of solving 

(3.34) W' = (A(x, A) - fis+)W, 

from X = +00 to X = for W"^ := A • • • A , where the hnear operator A is defined 
by the Leibnitz formula 

(3.35) AWi A • • • A VFfc := {AWi A W2 A • • • A W^) + • • • + (PFi A • • • A AWk) 

with eigenvectors and eigenvalues TZ = ri A • • • A r^, fi = ai + • • • + a^, where rj and aj 
are eigenvectors of A, and fis+ is the sum of the eigenvalues of A^ associated with the 
/c-dimensional stable subspace 5+, to obtain a solution asymptotic to an eigenvector TZs+ 
of A+ obtained as the wedge product of a basis of 5+; solving the symmetric equation from 
x = — cxDtox = for a solution asymptotic at —00 to an eigenvector TZs_ of A_ associated 
with the unstable subspace S- of A^; then evaluating the Evans function following (jl.6p as 

(3.36) Z)(A) := ((W+ AW-)U.=o), 

W~ := Wi A • • ■ /\W~_f^, where (•) denotes coordinatization in the standard Euclidean basis, 
i.e. 7] =: {rj){ei A ■ ■ ■ A e„) for an n-form r/H 

Equation (j3.34p is of the form (j3.ip considered in Section 13.11 so that we may apply 
our just-developed theory. Moreover, as ^13^ and respectively are the smallest real 
part and largest real part eigenvalues of A+ and A_, we are in the more favorable one- 
sided, or "Dirichlet" case S-t = C^* suitable for shooting methods, where A'^-i- := (^) and 

:= („"^) are the dimensions of the systems for and W~. We may thus approximate 
the solution as described in Section [3.21 bv solving a pair of finite difference schemes (j3.12p . 
on [0, M] and [0, — M], respectively, discretized asO = xo < xi < ■■■ < xj = M and 
= xo > X-i > ■ ■ ■ > x_j = — Af, with Dirichlet boundary conditions 

(3.37) wj = ns+, WZj = ns_, 
determining thereby a numerically approximated Evans function 

(3.38) V^'^\X):={iW+ AW,)), 

where h := (/i_j, . . . ,hj) is the vector of mesh blocks used to discretize [— M, M]. 

In the above discussion, we have implicitly assumed the consistent splitting hypothesis 
of |AG J| : that the dimensions of the stable subspace of ^+(A) and the unstable subspace of 
A- (A) sum to full rank n on the subset of A under consideration. By standard considerations 
|Hel \GZ\ IZlj . the "region of consistent splitting" on which this holds typically includes 
the component of real +00 in the complement of the essential spectrum of the associated 
linearized differential operator L of (jl.ip whose point spectra the Evans function is designed 



^See Section |6]3]T] for numerical prescriptions of TZs±- 
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to determine. However, as pointed out in |GZl IZH] . it is sometimes useful to extend this 
region of investigation and study also eigenvalues embedded in the essential spectrum. 

Following jGZl IZHj . we thus consider the problem in a slightly more general setting, 
substituting in place of consistent splitting the following gap assumption. 

Assumption 3.4. On A C C, the spaces 5"+ and S- are invariant subspaces of and 
A-, analytic in A, with dimensions k and (n — k). Moreover, the spectral gaps and u- 
defined as the maximum of the difference between the smallest real part of the eigenvalues 
of not associated with 5"+ and the largest real part of those associated with 5*+ and 
the maximum of the difference between the smallest real part of the eigenvalues of 
associated with S- and the largest real part of those not associated with ^.i., satisfy 

(3.39) > -e, j = ±, 

where > as in (|3.2p is the exponential rate of convergence of A to A± as x ±oo. 
Applying Theorem 13.91 in this context, we obtain the follow convergence result. 

Corollary 3.13. Under Assumptions lEM [13 and\3^ for fixed < 6 < 6, for M > 

sufficiently large and supj \hj\e^~^~"^^^ sufficiently small, where —9 < D < is less than 
m.miy±, there exist unique solutions of (j3.ip . (|3.37p determining an approximate Evans 
function T>^^'^{X) satisfying 

(3.40) \V^^^ -D\< C(e-^'*^ + sup |tj|) 

i 

uniformly on compact subsets of A, where D is constructed following (I3.36P from the solu- 
tions of (|3.34p guaranteed by Lemma \3.1l 9 is as in (|3.2p . and tj is truncation error. 

Proof. Immediate, observing that Assumption 13.41 implies Assumption 13.11 □ 

3.5.1 Mesh requirements, and computations in the essential spectrum 

So long as we remain in the region of consistent splitting (see discussion above Assumption 
13. 4p . i.e., away from the essential spectrum of the operator L whose point spectra we seek 
to study, we have 3?ctA > 0, with a simple eigenvalue at zero, from which we may obtain 
the stability property (j3.17p of (ii) needed for the convergence proof, with value /i = 0, 
even though the statement of (ii) is not strictly satisfied. Indeed, this situation holds for 
general difference schemes in case = C" whenever ^aA^ > and zero is a semisimple 
eigenvalue of A^, yielding the same convergence results stated in Theorem 13.41 and Corollary 
13.131 Moreover, the spectral gaps of Assumption 13.41 are identically zero, and D (by 
semisimplicity) may be taken as zero as well. 

Together with Remark 13.31 this shows that, for shooting methods based on an A-stable 
Cauchy solver (e.g., RK45), there is no requirement on the mesh size \hj\ beyond the 
requirement supj \hj\e^~^~'^^^^ < supj \hj\e^^^^ sufficiently small stated explicitly in the 



22 



Theorem (Corollary). This agrees with observations in \Br\ \BtZ\ IHuZH [BHRZI iHLZt ICHNZj 
in which centered exterior-product computations away from the essential spectrum are 
observed to require an extremely sparse mesh. 

On the other hand, for A inside the essential spectrum, one or more of the gaps 
becomes negative and ^aA± are no longer of one sign. As suggested by the simple com- 
putations of Remark 13.81 this might result in a much stricter requirement on | hj \ in order 
to satisfy condition (ii) (with /2 now strictly positive). Whether this is a real effect or just 
an artifact of our analysis is not clear, but this would be an interesting issue for further 
investigationlf] A second interesting question, for general centered schemes, would be the 
relation between the dimension of the kernel of A± and the size of the coefficient C in 
convergence estimate (I3.27P : we conjecture that they are roughly proportional, information 
that could be useful in comparing schemes. 

4 Stability of continuous orthogonalization 

We next address stability of the continuous orthogonalization method. Consider a solution 
O of the continuous orthogonalization system (jl.8p (i) restricted to the Stiefel manifold 

5 = {rj : 0,*^} = Ik}, such that 0, — > r2+ as x — > +oo. Evidently, the columns of 17+ are 
an orthonormal basis for an invariant subspace of = lim^^^+oo A{x). Of particular 
interest is the case arising in Evans function computations that is the stable or (at 
certain boundary points) a neutrally-stable subspace of A^. Linearizing (jl.Sp fi) about Q, 
we obtain the linearized system 

n' = {i- M*)An - nn*An - nn*An, 

for which the limiting constant-coefficient equation as a; — > +oo is 

(4.1) n' = m := (/ - n+n\)An - nnxAn+ - n+n*An+. 

We wish to assess the backward stability of ()4.ip with respect to general perturbations, 
both along the tangent manifold of the Stiefel manifold S and in transverse directions. 

Remark 4.1. Properly speaking, (j4.ip is linear in the pair of variables (17,0*) and not Q 
alone, since matrix adjoint is not a linear operation. Along the tangent manifold, however, 
it is linear as we shall see in 0, alone. 

4.1 Asymptotic stability in tangential directions 

The tangent manifold to S at consists of directions such that D(fl*i})Q_^_fl = 0, or 

(4.2) nin + n*n+ = o, 

®In any CEise, this would presumably be confined to shooting methods. 
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that is, for which fi^fi is skew-symmetric. As the Stiefel manifold is invariant under (jl.Sp fi). 
this is evidently invariant under ()4.ip . as direct computation readily verifies. It is spanned 
by the direct sum of eigenvectors ^l^K, K skew, in the kernel of C and eigenvectors 

(4.3) := (I - n+^\)r,al 

in the kernel of il^j., where Vj run through the eigenvectors of A transverse to and 
CTfc € C'^ are left eigenvectors of a defined by aQ.j^ := A^^-Q^^-. The former correspond to 
rotations within the same invariant subspace, the latter to perturbations outside S+. 
Under (j4.2p . we readily find that (j4.ip simplifies to 

(4.4) n' = m:={i-n+ni){An-na), 

which, for eigenvectors (14. 3|) yields 

Cfljk = ^jk{aj — ttk), 

where aj and ak are the eigenvalue associated with rj and cjfc- Thus, the eigenvalues along 
the Stiefel manifold are zero {k{k + l)/2-fold) and Uj — ((n — A;) x A: in total), where 
belong to and Oj to the compementary invariant subspace of Aj^. 

In particular, when S4. is the stable subspace of j4+, the Stiefel eigenvalues of C have 
either zero or positive real part, and so ()4.ip restricted to the tangent space, as claimed 
in the introduction, is (neutrally) stable in backward x, at least for the limiting system as 
X — > +00. 

Remark 4.2. Reduced equation (14. 4p is linear in Vt alone, since it does not involve f]*. 
4.2 Asymptotic stability in transverse directions 

Off the Stiefel manifold, we face the difficulty pointed out in Remark 14.11 that C strictly 
speaking is a linear function of i7 and il* , so that (|4.ip actually represents a pair of coupled 
equations, complicating calculations. However, we can sidestep much of this difficulty by 
noting that the remaining modes not already treated are of form where [3 S ^^^^ by 

a brief calculation satisfies 

(4.5) /j' = _(/3 + /3*)«, 

^+$7+ =: r2_|_a. Introducing TZ := K/3 := (l/2)(/3 + /?*), we thus have the linear system 



(4.6) 

which has block-triangular form 
(4.7) 



7^' = -TZa - a*n, 
13' = -27^a, 



n\' (M o\ fn 



(3 I \M Ol \p 
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where Mil := -TZa - a*n and MU := -ITZa. 

This has the /c^-dimensional kernel TZ = already identified in Section 14.11 and k'^ 
eigenvectors 

with eigenvalues ■= —{o-j + Ofc)) where Ij are left eigenvalues of a and aj the associated 
eigenvalues, which are exactly the eigenvalues of restricted to S+. 

In particular, when is the stable subspace of A^, the transverse eigenvalues of C have 
either zero or positive real part, and so the Stiefel manifold as indicated in the introduction, 
is (neutrally) stable in backward x, at least for the limiting system as x ^ +00. 

Remark 4.3. Alternatively, we may follow the simpler argument of the introduction to 
reach the same conclusion without finding explicitly the eigenmodes of the system. 



4.3 Convergence of the polar coordinate method 

Consider now the polar coordinate method, consisting of approximation of (11.81) on [— M, 0] 
and [0,M] by forward (resp. backward) Cauchy solvers, under Dirichlet boundary (i.e., 
Cauchy) conditions 

(4.9) W+ = (f^+,log7+), >V7 = (l^+,log7+), 
where are the numerical approximations of {Q^ ,logj±){xj). 

Corollary 4.4. Under Assumptions [13 \EM and for fixed < 9 < 9, for M > 

sufficiently large and supj \hj\e^~^~'^^^^ sufficiently small, where —6 < 9 < is less 
than min u± , there exist unique solutions of the discretization of (jl.Sp with boundary 
conditions (j4.9p determining an approximate Evans function D^'^'^{X) satisfying 

(4.10) \V^^''' -D\< C{e"^^ + sup |r,-|) 

j 

uniformly on compact subsets of A, where D is constructed following (j3.36p from the solu- 
tions of (j3.34p guaranteed by Lemma \3.1[ 9 is as in (j3.2p . and tj is truncation error. 

Proof. Equivalently, we must establish the analog of Theorem 13.41 This follows in routine 
fashion by decomposing the error equation for the numerical difference scheme into its linear 
and nonlinear parts and treating the nonlinear part along with the conjugation error from 
transformation into ^-coordinates together as source terms in a fixed-point equation in a 
combination of the discrete linear argument of Corollary 13.131 and the continuous argument 
of Lemma 12. H to obtain a contraction in the weighted ^""(Z"*") space defined by norm 
||VV|| := supj |Wje^^j |, similarly as in the standard proof of the Stable Manifold Theorem. 

We omit the straightforward but tedious details, except to mention one subtle point 
that will recur in later applications. Namely, the righthand side of (jl.Sp fi). considered 
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as a function on the complex- valued matrix Q, is not C^. For, it involves the operation 
of matrix adjoint, which in turn involves complex conjugation, a non-analytic function on 
complex arguments. Thus, we cannot immediately apply the above-described argument to 
(jl.Sp as written as a complex-valued ODE, but must instead first decompose it into real 
and imaginary parts, or, as we prefer to do, consider the doubled system 

= (i- nn)An, 

(4.11) . I 'J 

n' = nA*{i -QQ) 

in the pair of variables (17,17), with := fi*. With this (purely internal) change, the 
argument goes through as described to yield the claimed result □ 



5 Boundary-value algorithms 

Finally, on a more speculative note, we develop further some ideas of [Si IHuZlj regarding 
implementation of boundary-valued based Evans solvers for use in extremely large-scale 
systems, in the light of our new results. 



5.1 Sandstede's method 

We first describe (perhaps an imperfect translation of) Sandstede's original idea based 
on established projective boundary- value methods [Bel] . This consists, loosely speaking, 
of numerically solving the original, uncentered Evans system ()2.ip on [0, M] with mixed 
projective boundary conditions 

(5.1) U+Wf = 0, UoWo = am, m = l,...,k, 

where 11+ is the rank-(n — k) unstable eigenprojection of A-^-, Hq is a randomly chosen rank- 
k projection, and am, m = 1, . . . ,k are k random phase conditions determining a basis of 
k independent solutions VV". Here, as usual, denotes the numerical approximation of 
W"^{xj). For generic choices of Ho, a, this will give a numerically well-conditioned problem, 
so the procedure is to randomly select candidate values, then change these if the method 
does not converge after appropriate time. 

The solution of the boundary- value scheme for a given parameter value is then obtained 
by Newton iteration combined with continuation/path- following. Once the method is run- 
ning, initial guesses for Ho, a, and the solution itself may be chosen strategically based on 
the solution for nearby parameters, to improve conditioning/speed of convergence. 

The disadvantages of this method are two: (i) decay at plus infinity means we don't have 
control of the asymptotic behavior of solutions at +oo, making it difficult to impose the 
desirable property of analyticity in A; indeed, we prescribe solutions by phase conditions at 
X = 0, where we have no direct knowledge of the link to asymptotic behavior, (ii) (related) 

^ Note that applying a standard numerical difference scheme to the doubled system (|4.11|l yields an 
algorithm identical to what would be obtained by applying it to the original equation (|1.8|l (iV 
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prescription of random phase conditions at the origin is logistically compHcated, requiring 
additional error control/programming beyond just solution of the Evans ODE. 

Advantages of the method are the existence of a well-developed theory of conver- 
gence/error estimation for methods of this form, and a hoped-for dimensional advantage of 
iterative methods vs. shooting in the treatment of large systems. 

5.2 A polar coordinate-based method 

Following a suggestion of [HuZlj . we propose an alternative boundary- value scheme based 
on the polar coordinate method, using a Newton-based iterative boundary- value solveiH to 
approximate the solutions ,^~^) and {Q~ ,^~) of the doubled equations 

n' = {i- nn)Acn, n' = qa* {i - nn) , log 7 = Trace{nAcn) - Trace{nAcn)± 

on [0, zizM] and [— M, 0] with Dirichlet boundary conditions 

(5.2) {n,n,j){±M) = {n±,ni,^±), 

starting with the exact solution = (i^±,f^±,7±) at c = and continuing via the 

homotopy 

Ac:=A± + c{A-A±), c € [0, 1] 

to the desired solution of the full problem Ac = A at c = 1. 

This approach eliminates disadvantages (i)-(ii) of the standard approach and appears 
straightforward to code. Moreover, Corollarv 14.41 gives a rigorous convergence result, indi- 
cating at least theoretical feasibility. What remains to be seen is whether it is practically 
useful on the scale of interest, and how its performance compares with standard uncentered 
schemes as described in Section 15.11 We hope to address these questions in future work. 

Remark 5.1. As noted already in the proof of Corollary \4.4\ the use of doubled coordinates 
is necessary in order that the ODE be as a function of its arguments, since the matrix 
adjoint operation, since not analytic as a complex-valued function, is not . First-order 
smoothness is needed to apply Newton iteration, of which the first step is linearization. 

5.3 A conjugation-based method 

A possible drawback of the polar coordinate in numerically sensitive situations is its nonlin- 
earity. An alternative, still more speculative, linear solution would be to use a Newton-based 
iterative solver to approximate on [0, M] and [—M, 0] solutions and P~ of the conjuga- 
tion equations 

p' = AP := A±P - PA 
of (I2.13p . using projective boundary conditions 

n±(P±-/)(±M) = 0, UoP± = a 
*For example, MATLAB's BVP5P. 
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starting with initial guess =1 and using a similar homotopy 

Ac:= A± + c{A- A±), cg[0, 1] 
from A to its constant-coefficient limits A±, defining an Evans approximant simply as 
(5.3) P'^'*^(A) :=det(P+i?+,p-i2-)U=o, 

where are matrices whose columns are bases of the stable (unstable) subspace of A±^ 
Again, we have a rigorous convergence result, this time in the form of Theorem 13. 4| 
but it is not clear whether the scheme is practically useful, or if so how its performance 
compares to that of the previously mentioned schemes. Moreover, besides sharing difficulty 
(ii) of the standard method, it has the additional difficulty that the dimension of 114. may 
change with different A, necessitating still further modifications to the boundary conditions. 



6 Postscript: initialization of eigenbases at infinity 

For completeness, we describe, following |BrZl IHSZl IHuZH IZ2] . a simple and effective 
method for computing analytically-chosen initializing eigenbases at plus and minus spatial 
infinity. Combined with the integration methods described in the rest of the paper, this 
gives a basic working Evans solver that performs quite well in practice0 

6.1 Kato's ODE 

Denote by 11+ and n_ the eigenprojections of yl+ onto its stable subspace and onto 
its unstable subspace, with A± defined as in (jl.2p . Assume as in the introduction that the 
dimensions of the stable and unstable subspaces are constants k and n — A; on the desired 
region of investigation A E A and sum to n (the "consistent splitting condition" of |AGJj ). 
By standard matrix perturbation theory, II-i- are analytic in A for A G A [K]. Introduce the 
complex ODE 

(6.1) H = n'i?, R{K) = R^, 

where ' denotes d/dX, A* € A is fixed, 11 = II-i-, and R = R± with and R n x k and 
n X (n — k) complex matrices, and R^, is full rank and satisfies n(A*)i?* = i?*: that is, its 
columns are a basis for the stable (resp. unstable) subspace of A^ (resp. A_). 

Lemma 6.1 ([K |IZ2j ). There exists a global analytic solution R of (j6.ip on A such that (i) 
ranki? = ranki?*, (ii) UR = R, and (Hi) UR' = 0. 

Proof. As a linear ODE with analytic coefficients, (16. ip possesses an analytic solution in a 
neighborhood of A*, that may be extended globally along any curve, whence, by the principle 

^Described further in Section [S] 
Optimized versions may be found in the STABLAB package developed by J. Humpherys. 
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of analytic continuation, it possesses a global analytic solution on any simply connected 
domain containing A^, [K]. Property (i) follows likewise by the fact that R satisfies a linear 
ODE. Differentiating the identity II^ = n following [Kj yields nil' + n'll = 11', whence, 
multiplying on the right by 11, we find the key property 



From (j6.2p . we obtain 

{UR - R)' = {U'R + UR' - R') = H'R + (H - 1)11' R = liH' R, 
which, by nil'II = and = 11 gives 

(Hi? - R)' = -nn'(nfl - r), {hr - r){k) = o, 

from which (ii) follows by uniqueness of solutions of linear ODE. Expanding UR' = niT'i? 
and using IIR = R and HH'n = 0, we obtain UR' = lili'IiR = 0, verifying (iii). □ 

Remark 6.2. Property (in) indicates that the Kato basis is an optimal choice in the sense 
that it involves minimal variation in R. It is also useful as a direct characterization of the 
Kato basis independent of (j6.ip ; see IHS^ \BDG^ or Example \6.6\ below. 

6.2 Numerical implementation 

Choose a set of mesh points Xj, j = 0, . . . , J along a path F C A and denote by Hj := n(Aj) 
and Rj the approximation of R{Xj). Typically, Aq = Aj, i.e., F is a closed contour. 

6.2.1 Computing Ilj 

Given a matrix A, one may efficiently (~ 32n^ operations; see |GvLl ISBj ) compute by 
"ordered" Schur decomp ositiorl^. i.e., Schur decomposition A = QUQ^^, Q orthogonal 
and U upper triangular, for which also the diagonal entries of U are ordered in increasing 
real part, an orthonormal basis 



of its unstable subspace, where Qk+i, ■ ■ ■ , Qn are the last n — k columns of Q, n — k the 
dimension of the unstable subspace. Performing the same procedure for —A, A*, and 
—A* we obtain orthonormal bases Rg, L^, Lg also for the stable subspace of A and the 
unstable and stable subspaces of A*, from which we may compute the stable and unstable 
eigenprojections in straightforward and numerically well-conditioned manner via 



(6.2) 



nn'n = o. 



R. 



■u 



{Qk+1, . . . ) Qn) 



(6.3) n, := Rs{LIRs)-'l:, U, 

Applying this to matrices A^ := ^-i-(Aj), we obta 
we consider Hf as known quantities. 
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Supported, for example, in MATLAB and LAPACK. 
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Remark 6.3. It is tempting to instead simply call an eigenvalue-eigenvector solver and 

r 1* 

express II^ = ^ where rj, Ij are left and right eigenvalues associated with stable eigen- 
vectors. However, this becomes ill-conditioned near points where stable eigenvalues collide, 
as frequently happens for cases resulting from other than simple scalar equations. 

6.2.2 First-order integration scheme 

Approximating n'(Aj) to first order by tiie finite difference (Hj+i — Hj) / {Xjj^i — Xj) and 
substituting this into a first-order Euler sclieme gives 

Rj+i = Rj + (Aj+i - Xj)-^ r'^ji 

or Rj+i = Rj + Hj+iRj — Iljiij, yielding by tlie property li-jRj = Rj (preserved exactly by 
the scheme) the simple greedy algorithm 

(6.4) Rj+i = Uj+iRj. 

It is a remarkable fact [Z2\ (a consequence of Lemma 16. ip that, up to numerical error, 
evolution of ()6.4p about a closed loop Aq = Aj yields the original value Rj = Rq. 

6.2.3 Second-order scheme 

To obtain a second-order discretization of (|6.ip . we approximate Rj+i—Rj AXjIl'j^-^^^2^j+i/2i 
good to second order, where AAj := Aj+i — Aj. Noting that Rj_^_i/2 ~ Ilj-^-i/2Rj to second 
order, by ()6.4p . and approximating Ilj_^_i/2 ~ ^{Hj+i + nj), also good to second order, and 
~ i^j+i ~ nj)/AAj, we obtain, combining and rearranging, 

Rj+i = Rj + ^(n,+i - n,)(nj+i + n,)i?,. 

Stabilizing by following with a projection Hj+i, we obtain after some rearrangement the 
reduced second-order explicit scheme 

(6.5) Rj+i = n,+i[/ + ^Uj{i - n,+i)]i?,. 

This is the version we recommend for serious computations. For individual numerical ex- 
periments the simpler greedy algorithm (j6.4p will often suffice (see discussion, [Z2]). 

Remark 6.4. Arbitrarily higher-order schemes may be obtained by Richardson extrapolation 
starting from scheme ()6.4p or (j6.5p ." see fZ^. In practice, this does not seem useful. 

6.3 Initialization of Evans function ODE 

Finally, we describe the conversion of analytic bases -R±(A) into initial data for the centered 
exterior product or polar coordinate method. 
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6.3.1 Centered exterior product scheme 

Denote = {Rf , ■ ■ ■ , R'l) and R^ = , . . . , Then, the initiahzing wedge prod- 

ucts TZs± of Section 13.51 are given simply by 

(6.6) ns+:=Rt ^ Rl and 7^s_ ■.= R^ R-_^. 

6.3.2 Polar coordinate scheme 

For each A E A, we may efficiently compute matrices r2-|-(A) whose columns form orthonor- 
mal bases for S± , by the same ordered Schur decomposition used in the computation of Il± . 
This need not even be continuous with respect to A. Equating 

n+a+{X) = R+{X), 17_a_(A) = R-{X), 

for some a±, we obtain 

a+(A) = rj;i2+(A), a_(A) = n*_R-{X), 

and therefore the exterior product of the columns of R± is equal to the exterior product of 
the columns of times 

(6.7) 7±(A) :=det(l^*i?)±(A)). 
Thus, we may initialize the polar coordinate ODE (jl.Sp with 

(6.8) n = n±, 7 = det{n*R)±. 

Remark 6.5. The wedge products represented by polar coordinates {'j,Q)±{X), with ^± 
defined as in (j6.7p . are the same products TZs± defined in ()6.6p . In particular, they are 
analytic with respect to X, though coordinates 7 and $7 in general are not. 

6.4 Error control 

As the integration of Kato's ODE is carried out on a bounded closed curve, standard error 
estimates apply and convergence is essentially automatic, and we shall not discuss it. In 
applications, we are often interested in determining the winding number of D about such a 
curve. For this purpose, following [Bt\ IBrZ] . we introduce a simple a posteriori "Rouche" 
check limiting the step-size in A by the requirement that the relative change in -D(A) be less 
than a conservative 0.1. (By Rouche's Theorem, relative error less than one is sufficient to 
obtain the correct winding number.) In contrast to integration in x of the Evans system, 
integration in A of the Kato system is a one-time cost, so not a rate-determining factor in the 
performance of the overall code. However, the computation time is sensitive (proportional) 
to the number of mesh points in A, so this should be held down as much as possible. 
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6.5 Finer points: two exceptional cases 



We conclude by pointing out two commonly occurring cases that can give trouble if not 
expected, and describe some practical resolutions. 



6.5.1 Behavior near the origin 

For the linearized operators L arising in the study of stability of traveling-waves of certain 
systems such as viscous conservation laws or Cahn-Hilliard and nonlinear Schrodinger equa- 
tions, the point A = is embedded in the essential spectrum of L. In computing a winding 
number around some bounded portion of the set {5?A > 0} of possible unstable eigenvalues, 
we must pass through or near this value, at which the spectral gap (see discussion above 
Assumption 13. 4p between stable and unstable subspaces of A± goes to zero. In this case, 
the eigenprojections II-i- lose their characterization as stable (resp. unstable) eigenprojec- 
tions of A±, so must be computed in a different way than the ordered Schur decomposition 
described above. Worse, they may lose analyticity, possessing a branch singularity, at A = 0. 

To avoid the former problem, we typically just compute near but not at A = 0. However, 
this leads to occasional uncertainty/bad results near the origin and should probably be 
improved in the analytic case by instead computing II-i- at points within or on the essential 
spectrum boundary of L by analytic extrapolation from values at points outside. This 
is an important practical area for further algorithm developmental The latter problem, 
concerning behavior near a branch singularity, is discussed in the next subsection. 



6.5.2 Behavior near a branch singularity 

For certain problems, especially those involving additional parameters, e.g. multi-d |HLyZ2| 
or families of one-dimensional waves that pass through characteristic values |BHZj . there 
may appear for certain parameters branch points in the eigenvalues of A± as a function of 
A, at which II-i- therefore blow up [K]. This requires some adustment in order to restore 
good behavior. 

Example 6.6. A model for this situation is the eigenvalue equation for a scalar convected 
heat equation Xu + r]u' = u" with convection coefficient r] passing through zero. The coeffi- 
cient matrix for the associated first- order system is 



(6.9) A :-- 



'0 r 

Then, the stable eigenvector of A determined by Kato 's ODE (j6.ip is 
(6.10) Riv, A) := ^^y^l^y/, (l, -v/2 - vWTA ) , 

which, apart from the divergent factor |^2^4^^ji/4 ; is a smooth function of 7/^/4 + A. 

^^This would also allow computations within the essential spectrum, up to now not systematically carried 
out (though see Br for some preliminary results in this direction). 
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Proof. By straightforward computation, fJ.±{X) := ^{rj/l+y^rj'^/i + A and V± := {1, ij.±{X))'^ 
are eigenvalues and eigenvectors of the matrix A of (j6.9p in Example 16.61 The associated 
Kato eigenvectors are determined uniquely, up to a constant factor independent of A, 
by the property that there exist corresponding left eigenvectors such that 

(6.11) (y • y)^ = constant, {V-V)'^ = 0, 

where " ' " denotes d/dX; see Lemma |6. If iii). 

Computing dual eigenvectors = {X+ii^)~^{X, fi±) satisfying (V- V)^ = 1, and setting 
= c±V^, = V^/c±, we find after a brief calculation that (jO.lip is equivalent to the 

complex ODE 

which may be solved by exponentiation, yielding the general solution c±(A) = C(r/^/4 + 
^-j-i/4_ Initializing without loss of generality at c±(l) = 1, we obtain (jG.lOp . □ 

Remark 6.7. It is straightforward to generalize by the same method the computation of 
Example \6.6\ to branch singulariteis of general order s. 

The computation of Example (j6.6p indicates that the Kato basis blows up at A = as 
(ry^ + 4A) as rj crosses the characteristic value r/ = 0, hence does not give a choice that 
is continuous across the entire range of parameters. However, the same example shows that 
there is a different choice (1,— t//2 — y'rf/A + X)'^ that is continuous, possessing only a 
square-root singularity. We can effectively exchange one for another, by rescaling the Kato 
basis by factor (r/^ + 4:X)^^^. See |BHZ| for examples/further details. 

This issue is mainly important in problems with parameters, or possessing branch singu- 
larities, but can also arise for a problem without parameter or singularity that happens to 
lie near a related problem with branch singularities. In such a case the "invisible" branch 
singularity could serve as an organizing center directing the Kato flow without the user 
being aware of it. Thus, it is important to be alert to this possibility. 
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